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ABSTRACT 

We derive a semi-analytical extension of the spherical collapse model of structure formation that takes 
account of the effects of deviations from spherical symmetry and shell crossing which are important in 
the non-linear regime. Our model is designed so that it predicts a relation between the peculiar velocity 
and density contrast that agrees with the results of N-body simulations in the region where such a 
comparison can sensibly be made. Prior to turnaround, when the unmodified spherical collapse model 
is expect to be a good approximation, the predictions of the two models coincide almost exactly. The 
effects of a late time dominating dark energy component are also taken into account. The improved 
spherical collapse model is a useful tool when one requires a good approximation not just to the 
evolution of the density contrast but also its trajectory. Moreover, the analytical fitting formulae 
presented is simple enough to be used anywhere where the standard spherical collapse might be used 
but with the advantage that it includes a realistic model of the effects of virialisation. 

Subject headings: Cosmology: Theory, miscellaneous. Relativity. Galaxies: general. 



1. INTRODUCTION 

The S pherical Collap se Mod el (hereafter SCM) devel- 
oped by iGunn fc Gottj (|1972f ) is perhaps the simplest 
model for the evolution of non-linear structure, and yet 
it has been shown to be remarkably successful when cor- 
rectly interpreted. However, despite the SCM's many 
success, it is ultimately flawed since it predicts that any 
overdensity of matter collapses to a singularity in a finite 
time. Additionally, making the assumption of spheri- 
cal symmetry, whilst simplifying, means abandoning the 
many interesting and important aspects of structure for- 
mation that result from deviations from spherical sym- 
metry. Indeed, these deviations play a crucial role in ul- 
timately halting the collapse of the overdensity and the 
formation of virialized structures. 

The usual approach to virialisation in the SCM is to 
put it in by hand] the collapse is simply halted once the 
virial radius has been reached. This procedure, within 
an Einstein-de Sitter Universe, leads to the result that 
bound structures are formed when the non-linear over- 
density is about 178 or, equivalently, the linear overden- 
sity is approximately 1.68 ^. Despite the fact that these 
figures are not too far away from what is actually ob- 
served in N-body simulations, the ad hoc nature of this 
approach means the SCM cannot be used to predict the 
precise manner in which the overdensity evolves; more- 
over, the SCM's prediction for the peculiar velocity be- 
c omes virtually useless shortly after turnaround. 

[Engineer et al.l 1)20001 ) proposed a different and better 
motivated way in which the SCM could be extended to 
include virialisation. Their idea was to alter the standard 
SCM evolution equation for the mean density contrast, 
(5, in such a way that stable structures would form. Their 
modified evolution equation was constructed by adding 

^ In dark energy scenarios these v alues change slightly depend- 
ing o n the model (Percival 2005; Deb nath et al.ir2006t iBagla et al.l 
[20051: IWeinberg fc Kamionkowskiii200f i 



a Taylor series in 1/5 to the standard equation: the idea 
being that these additional terms would encode all the 
effects due to shell-crossing and deviations from spheri- 
cal symmetry that occur for large 5. The coefficients of 
the terms in Taylor series were chosen so as to provide 
a good approximation to the statistical density-contrast 
found from N-body simulations. Whilst their improved 
SCM was found to agree fairly well with data from N- 
body simulations for 5 > 15, their solution for 5 is in- 
accurate in the linear regime, where spherical symmetry 
and hence also the unmodified SCM are expected to be 
good approximations. 

In this art i cle w e take a similar approach to 
[Engineer et al.l ()200G[ ) and add terms to the standard 
SCM evolution equation so that for large 5 the evolu- 
tion of the density contrast and the peculiar velocity 
agrees with the data from N-body simulations. Impor- 
tantly though, the analytical solutions, which we find, 
remain valid as (5 — > and agree almost exactly with the 
unmodified SCM prior to turnaround. 

2. MODELLING NON-LINEAR STRUCTURES 

We consider a spherical overdensity embedded in 
a background universe described by a Friedman- 
Robertson- Walker (FRW) metric. In a spherically- 
symmetric system and in the absence of shell-crossing, 
the mass, M, inside each comoving spherical shell, with 
physical radius R{t), remains constant. The mean den- 
sity contrast inside shell therefore scales as: 



l + <5 = A 



a' it) 
R{tf 



(1) 



where a{t) is the scale factor of the background FRW uni- 
verse, and A is constant on each shell. Eq. JT]) can also 
be used to define R{t) for shells of constant M when de- 
viations from spherical symmetry occur. In these cases, 
howeve r, the physical meanin g of R{t) is less clear. Fol- 
lowing [Engineer et al.[ ()200Gf ) we continue to think of 
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R{t), as defined by Eq. ([T]), as the effective "radius" 
of each, shell of constant M. 
The peculiar velocity, hsci is defined by: 



i-sc 



1 dln(l + (5) 
3 din a 




(2) 



where H = a/a is the Hubble parameter. In a matter 
dominated Universe the unmodified SCM predicts hsc = 
hsci^)- Our key assumption in what follows is that this 
remains the case even when deviations from spherical 
symmetry are taken into account. This assumption is 
supported by numerical s tudies of structure for mation 
such as that conducted bv [Hamilton et aTl (jl991f ) which 
we discuss below. 

For the moment we take the background universe to 
be matter dominated, however in Section[5]we generalize 
our results to non-linear overdensities in more realistic 
universes which have recently transitioned to an epoch 
of dark energy domination. We define a new coordinate 
by: 



1 



2(l-cos?7)3' 

for some function T{ri). Now M is constant on each 
shell and so, taking t ~ ti, R = Ri and 6 ~ 5i &i some 
instance, we have 



M 



2Rl{l + 5.,) 2R^{1 + S) 



9t? 



9t2 



By considering this relation together with Eq. ^ we see 
that: 



R = 



t- 



R^{l + 5,) 

3U{l + 5, 



(1 - COs(7y)) 



where j{r],Si) is some function of r/ and Si. If hsc — 
hsciS) then j{ri,Ri) — iifi), and without loss of gener- 
ality we can set j{vi) = 1. This gives: 



sm77 



dr 



2 (1 -cos 77)1/2 dT' 
where T(jf) = i] — sinry, and: 



(4) 



i? = i?(^)^:^iii±il)(l_cos,?), (5) 
* = *i^) = -ly^ — ^W' ■ 



4,5: 



(6) 



We reduce to th e standard SCM in t he limit where T = 
T{r]) = rj — sin rj ()Padmanabhan|[2002f ) . and it is clear that 
hsc — hsc{S) in the unmodified SCM. It is important 
to stress that if hsc = hsci^) then Eqs. ([5]) and © are 
the most general solutions for R and t. 
The radial acceleration of each shell is found to be: 

d^R 



GM 



r'-2+sinr7(l-cosr7)|^) , (7) 



where T' = dT/dr, and T" = d^T/dr^. As should be 
expected, when T' = 1 we reduce to the standard SCM 
equation for R^tt- 



When deviations from spherical symmetry and the 
leading order effects of a gradie nt in the veloc it y dis - 
persion are taken into account, [Engineer et al.l ( 20001 ) 
showed that one should have: 



d^R 

dt2 ^ 



GM H^R 



i?2 

GM 
'~R^ 



(8) 



3(1 + 5) 



S 



where S{a,x) — a'^{a^ — 2il^) + f{a,x); and Q"^ re- 
spectively quantify the shear and rotation of the fluid; 
/(a, x) contains the lowest order contribution from veloc- 
ity dispersion terms. Importantly, no matter what form 
S{a, x) takes, if, as we have assumed hsc = hsc{S), then 
we must have S = S{S) and by comparing Eqs. dZl) and 
([8|) we can clearly see that: 



3(6) 



3(1 + ^) 



T'-^ - 1 + sinr;(l - cost?) 



rp/l 



In principle, the form of both S{6) and hence T{r]) can 
be found using the results of N-body simulations. Unfor- 
tunately, however, making the required comparison with 
simulations is not as straightforward as one might expect 
it to be. This is because the results of such simulations 
are given in terms of the statistical properties of the mat- 
ter distribution rather than in terms of the mean density 
contrast, S, and peculiar velocity, hsc- The statistical 
properties in question are the averaged two point corre- 
lation function, ^, and the averaged pair velocity, h{a, x), 
which are given as defined by: 



£^{x,a)x^ dx; h{a,x) = 



{v{a,x)) 



(9) 



where ^ is the two-point correlation function and is de- 
fined to be the Fourier transform of the power spectrum, 
P{k). The assumption that h{a,x) depends on a and x 
only through £ i.e. h(a,x) = h(^(a,x ) ) is commo n in the 
literature fsee [Engineer et al.l pOOOl ). iMoe et al.l (|1995D 1 
and it appear s to have been confirme d by numerical sim- 
ulations (see iHamilton et al.l (|199lD . iPeacock fc Doddj 
([1996)). 

The results of, for example, IHamilton et al.l (119911 ) 
can be used to cons truct the fitting formula for h{^) 
(|Engineer et al.|[2000l ). However, before we can make use 
of such a forrnula, we must relate the statistical quanti- 
ties £ and /t(£) to 6 an d hsc(S). It is well know n that 
( Padmanabha nI (ll996DlPadmanabhanl (2002) , Peebles 
ll980) .Padmanabhan fc Eng,irieej (|1996>).Popolo et al. 
(2001) jPadmanabhan fc Rayl (|2006l )). on scales smaller 
than the size of the collapsing objects and around high 
density peaks: 

p-Pfc(i + 0- (10) 

It follows that, in the non-linear regime, we have S ~ ^. 
This relationship be tween S and £ was also used by 
lEngineer et al.l ()2000f ). although it was, as it is here, the 
weakest part of the whole analysis. We only require that 
(5 « £ hold where it is expected to be a good approxi- 
m ation i .e. S 15. As (5 — > qo we assume that 5 ^ ^. 
iPeebled (|1980( ) showed that /i(£) satisfies: 



h = 



1 1 



3 1 + £dlna' 



(11) 



3 



thus if ^ « (5, then by comparing Eqs. ([2]) and (fTTj) . we 
see that h ~ hsc- It must be stressed that this second 
relation between hsc and h is only valid if ^ « 5, however 
whenever itjdoes hold it implies that hsc, is given by a 
function of S alone, i.e. hsc — hsci^). 

The assumption that ^ w ^ breaks down for small S. 
Fortunately, when 6 is small, the unmodified SCM pro- 
vides an accurate model. For times _i ^ ti, the unmodi- 
fied SCM predicts that hsc = hsc{S)- Our key assump- 
tion that hsc — hsc{^) is therefore expected to hold for 
almost all S. The assumption does break down for t ~ 
however this is entirely due to the decaying mode in 6 
which is negligible for t ^ ti. 

3. CONSTRUCTING AN IMPROVED SCM 

If hsc = hsc{^) then all properties of a modified SCM 
are encoded in a single function T{r]). Our aim is to com- 
bine the unmodified SCM and data from N-body simu- 
lations to find a fitting formula for T{'i]) that results in 
accurate predictions for hsci^) in all regimes. 

Prior to turnaround we expect the standard SCM to 
be accurate and hence T{r]) « 77 — sin 77. Furthermore, for 
the unmodified SCM to be accurate, at leading order, in 
the linear regime we must have T(77) ~ 77 — sinrj + 0(77'^) 
for small rj. When 6 > 15, we expect hsc{^) ~ ^(C) 
and ^ ~ S, and we may use h{^) to extract the large S 
form of T'(77). We describe how this is done below. In 
a matter dominated universe, the linearly extrapolated 
mean two-point correlation function, , scales as oc 
cx T''/'^, therefore as (5 — > 00: 



with T « Tsim = Vsim " smTjsim] Vsim is defined by 



(i + 0«(i + ^) 



T(77)2 



2 (1 -cos(77))3' 



and so 



(1 - cos 77) W (1 - COSIJsim) 



(12) 



(13) 



(n-e>/3' 

where A is a constant and we tre at it as a parameter 
to be fitted. [Hamilton et al.l ()1991l ) found the following 
fitting formula for ^im(^): 

- ^/ 1-H0.0158e^ + 0.000115C^ V'^ 

Ku-a - ^ 1^ 1 _,_ 0.926^2 - 0.0743^3 -f 0.0156^* ) ' ^ ' 

We define 7700 — limj_,g^ 77. By taking ^ cxd we find 
that: 



A = 2.2668(1 - cos 7700) = 2.2668 



Rta 



(15) 



where we have used cos 7^00 — ^^'^Rvir/ Rta which follows 
from Eq.(l5|); Rvir is the radius of the shell at virialisation, 
and Rta is its radius at turnaround. If spherical symme- 
try is assumed then the virial theorem in an Einstein 
de Sitter Universe gives Rmr = Rta/"^- This relation 
is generally used when virialisa tion is placed b y hand 
into the SCM. For comparsion, iHamilton et~aI1 (|199lD 
found that Rta/ Rvir ~ 1-8 from their simulations. We 
treat Rvir /Rta as a fitting parameter. Eq. (jl3p provides 
77 = 77(^ « (5), and T(77) may now be fomid using Eq. 

dni): 



1/2 



(16) 



We fit for the parameter A (or equivalently Rta/ Rvir) 
by considering some important physical constraints on 
the behaviour of r(77). Since the effect of deviations from 
spherical symmetry is to slow down the collapse of the 
overdensity, it follows from Eq.([7]) that: 



rpll 

r'^ - sin 77(1 - cos 77)— > 1. 



(17) 



Furthermore, the unaltered SCM should be a good ap- 
proximation up to around turnaround (when 77 = tt) and 
so we must have that T' Ri 1, T" < T' for 77 < tt. As 
T — > Voo — sin 7700 , we must have T — > 00 and so T' > 
as 77 — » 7700. Eq. I|17p implies that, for 77 > tt, we can- 
not have < T' < 1 and T" < 0, and so it follows 
that T' > 1 for 77 > TT. It follows that T > t every- 
where. This condition implies that we should choose 
Rmr/Rta, and hence A, so that Tgi^ is always greater 
than Tsim = rjsim - siin]sim for f > 15, which roughly 
corresponds to 77sim ^ 3.8, Tsim ^ 4.4. Moreover, since 
we want T ^ t as 77 — > 0, we choose Rmr/Rta so that, 
at the minimum of Tsim, Tsim = Tsim- This requirement 
gives: 

^ = 0.5896. 

Rta 

We now use Tsim to find a fitting formula for 2^(77) that 
agrees with the unmodified SCM at early times (i.e. T 
7] — sin 77 for small 77) . 

4. RESULTS 



We find that the fitting formula: 



T(77(r)) 



3.468(r/ -T)-i/2exp 



(1 + 0.8(r/ - r)i/2 _ o.4(t/ - t)) 



(18) 

where Tf — 5.516, provides an excellent fit to the form of 
T(n(T )) derived from the simulations of IHamilton et aP 
(|1991h . i.e. T « Ts^m, in the range 6 > 15. It also 
provides an evolution of S that matches up smoothly to 
the one predicted by the standard SCM in the region 
where we expect it to provide an accurate approximation, 
i.e. prior to turnaround (5 < 5. 

In FIG [1] we plot Tsim against Tsim , and our fitting 
formula for T{ri{T)) against r. We see that, in the range 
r>4.4 where we expect ^ « 5 and hence Tsim ~ 
T , the fit is indeed very good. The parameters in the 
formula for T{rj{T)) have been chosen so that as r — > r/, 
the leading order terms in the asymptotic expansions of 
both T(77(t)) and TsimiVsimir)) match. 

In FIG El we use Eqs. Q a nd ffl to plot hsc{5\ 
and the fitting formula found by IHamilton et al.l (|1991l ) 
to plot h{£). We also plot the unmodified SCM pre- 
diction for hsc{S). It is clear from this plot that our 
fitting formula for T gives an hsc{S) that is an excellent 
approximation to /i(^) (provided ^ « (5) in the region 
6 > 20. As ^, ^ — > 00, the curves hsc{S) and h{^) have 
the same leading order asymptotic behaviour. In the re- 
gion S > 15, our model gives an hsci^) that is always 
within 3% of the fitting for mula for h{^) derived from 
simulations (|Hamilton et al.|[l991 ). Additionally, our fit- 
ting formula for hsc{S) agrees almost exactly with the 
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Fig. 1. — Plot of T/t versus t. The solid red line is T/r = 
Taim/Taim against T = Taimy whereas the dashed black line shows 
T(t) /t i n our improved SCM and uses the fitting formula for T(t), 
Eq. II18I I. to evaluate T/t. As required, we see that in the region 
T > 4.4 <5 ^ 15, our improved SCM is a very good approximation 
to the simulation data. For small values of r, we see that T/t — > 1, 
as is required, in our improved SCM. 



Fig. 3. — Plot of how the non-linear mean density contrast, 5, 
depends on the linear one, <5;i„. The solid red line is the improved 
SCM developed here, the thick dashed green line is the standard 
SCM with virialisation put in by hand and the thin dashed black 
line is the standard SCM without virialisation. The dotted blue line 
shows the asymptotic behaviour of <5(<5;i„) that we would expect 
if Rta /Ryir = 1.8 as seen in the simulations of IHamilton et aU 
llT99lh . 




1+5, 1+^ 

Fig. 2. — Relationship between h and S in different models. 
The solid red line shows hsci^) in onr improved spherical col- 
lapse model. The dotted black line is hgci^) in th^ standard 
SCM, and the dashed blue line shows as seen in simulations 
IIHamilton et al.1ll991h . As required, our improved SCM gives an 
hgc that agrees with h{^) very well in the region where we expect 
5 ^ i.e. 5 > 15. We also note that hsci^)} in the improved SCM, 
has the same asymptotic behaviour, to leading order, as Rj S) 
in the limit 5 — > oo. Importantly, our improved SCM model also 
agrees with the standard SCM in the region where it is expected 
to provide a very good approximation i.e. prior to the epoch of 
turnaround, 5 < 5.6. 



predictions of the unmodified SCM prior to turnaround 
6 < 5.6. We found a best fit value of Ryir/Rta — 0.5896, 
wliicli is fairly close to the value of 0.5 that is gener- 



ally used in the spherical collapse model, despite the 
fact that we did not constrained it to be so. For com- 
parison, the improve d spherical collapse developed by 
[Engineer et"aLl (|2000D gave Rvir/Rta « 0.65, and N-body 
simulations have been found to support Rvir/Rta ~ 0.56 
(|Hamilton et al.|[l99lD . Our value of Rvir/Rta is there- 
fore a better approximation to the value found from sim- 
ulations than that generally used in the u nmodified SCM 
and th at found in the model developed bv lEngineer et ah! 
(|2000l) . 

In our improved SCM the linear density contrast is 
given by: 

When Siin = 1.6865, which corresponds to the instant of 
collapse in the unmodified SCM, we find 6 = 54.65 rather 
than the SCM's value of 178. We find that d = 200 cor- 
responds to Siin — 2.286. We compare the form of S{Siin) 
found in our model to that predicted by the unmodified 
SCM model, with virialisation put in by hand^ in FIG. 
m We also show the divergent behaviour of S{Siin) in 
the standard SCM without virialisation, and the asymp- 
totic behaviour of S(Sii n ) if Rta = l-8 Rvir as suggested 
by simulations (Hamil ton et al.lll991h . The two major 
advantages of our improved SCM over the unmodified 
version are clearly visible in this plot: 

• At late times the improved SCM provides a bet- 
ter approximation to the behaviour seen in N-body 
simulations than the unmodified SCM does. 

• S(Siin) is smooth in the improved SCM. This is not 
the case in the unmodified SCM when virialisation 
is put in by hand. The improved SCM therefore 
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provides not only a good approximation to the evo- 
lution of density contrast, S, but also to that of 
dS I At and hence to that of the peculiar velocity. 

5. INCLUDING DARK ENERGY 

We have constructed an improved, semi-analytical 
SCM whose prediction for the dependence of peculiar 
velocity on the mean density contrast concurs with that 
extrapolated from simulations in the non-linear regime. 
However, as mentioned above, we have done this for a 
matter dominated universe. It is well known, however, 
that the universe today is not matter dominated and a 
significant fraction of its total content is believed to be 
in the form of dark energy. Fortunately though, the ef- 
fects of this dark energy on the background universe are 
generally believed to only have become non-negligible rel- 
atively recently z < 1.8. This means the local evolution 
of overdensities with 6 > 100 today has been matter- 
dominated, to a very good approximation, right up until 
the epoch of matter-radiation equality. Therefore the 
evolution of the density, p, and radial velocity, of a 
shell of matter that is well into the non-linear regime to- 
day is, to a good approximation, the same in a matter 
dominated background as it is in a background where 
today Qm ~ 0.27 and Qde ~ 0.73 for which the equation 
of state parameter for the dark energy, w, satisfies the 
current astronomi cal bound oi w — — 1 ± 0.1 for z < 1 

less et ail I2006D . It is therefore a fairly straightfor- 
ward task to generalize our formulae for 5 and hsc to 
included background universes where Qm ~ 0.27 pro- 
vided that these are only applied to overdensities that 
are large today. We then find that: 



2(1 — cos 77)-^ 



(19) 



hscinrn,t) = i-a 



0.5 



' 1 + S dr , , 

— -snn— 20) 

2(1-0087;) 'dT ^ ' 



where t = rj — sinr/ as before, T{ri) is still given by Eq. 
dUl) and 

The equations for R{'q) and t{rf) remain unchanged. 

If dark energy has behaved similarly to a cosmological 
constant in the recent past {z < 1.8) then: 



/(«) 



r!,„ (sinh- 1 ( v/(l - a„)/f^™) ) 2 



-0.4 



Eqs. (|T9|) and (|20p combined with the fitting formula 
Eq. psp provide a very good approximation to the evo- 



lution of overdensities of matter in a realistic universe 
provided that 5 > 100 today. For smaller values of 5 
our results are only accurate for « 1. However for 
(5 < 15 the effects of deviations from spherical symmetry, 
which have been our primary concern in this article, are 
small enough to be ignored and the results derived using 
the unmodified SCM can be used with confidence. No- 
tice also that such fitting formulae and assumptions are 
strictly speaking only applicable to dark energy models 
that are not coupled to matter. If such a matter cou- 
pling is allowed then the whole process of structure for- 
mation, both in t he linear and non - hnear r egimes, may 
change (see e.g. iB rookfield et al.' ('2006'); iLaha v et al " 



1991D; 
12005); 



Mota k va n dc Bruck (2004) 

(2006); iMota fc Shawl 



Nunes & Mota 
|2007); Koivi sto fc Mota (2007a,jJ)) 



Maor & Lahav 
( 200a 



This said, in many 



dark energy models, the matter coupling is constrained 
by experiments to be very small (relative to the coupling 
between matter and gravity), and as a result any alter- 
ations to the process of structure formation are similarly 
small. 

6. CONCLUSIONS 

In this article we have extended the Spherical Collapse 
Model so that it takes account of the effects of devia- 
tions from spherical symmetry and shell crossing which 
are important in the non-linear regime. The key assump- 
tions that we^used when constructing oui^ model was that 
hsc = hsciS) and that for (5 > 15, (5 « ^. The latter as- 
sumption is probably the weakest link in the whole anal- 
ysis, although both assum ptions are found commonly 
in the literatur e (see Engineer et al.l (|2000( ): iMoe et all 
lfT995h : iPeeblej (fl98nl) V Our improved SCM predicts a 
form for hsc{S) that is consistent with the results of N- 
body simulations in the regime where a comparison can 
sensibly be made {S > 15) and with the unmodified SCM 
prior to turnaround. Analytical formulae for S and hsc 
in the improved model have been presented, and they es- 
sentially differ from the comparable formulae in the un- 
modified SCM only by the replacement of t = 77 — sin 77 
with T{ri), which is given by Eq. (jl8p . The improved 
SCM is therefore simple enough to be used anywhere 
where the unmodified SCM might be used but with the 
advantage that it includes a realistic model of the effects 
of virialisation. 
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